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We theoretically derive the amplitude equations for a self-propelled droplet driven by Marangoni 
flow. As advective flow driven by surface tension gradient is enhanced, the stationary state becomes 
unstable and the droplet starts to move. The velocity of the droplet is determined from a cubic 
nonlinear term in the amplitude equations. The obtained critical point and the characteristic velocity 
are well supported by numerical simulations. 
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I. INTRODUCTION 

Spontaneous motion or self-propulsion has been at- 
tracting attention in recent decades because of its poten- 
tial application to biological problems such as cell motil- 
ity [IHS] . These intensive studies have stemmed from the 
fact that mechanical properties of cells can be measured 
thanks to recent developments in visualization techniques 
[Hj . In addition, several model experiments showing spon- 
taneous motion have been carried out [7T-FT2] . These sys- 
tems consisted of relatively simple components such as oil 
droplets in the water (5]. Nevertheless, the droplets give 
the impression of being alive in that they move sponta- 
neously without being pushed or pulled, and they travel 
in straight lines, turn, and deform. 

Motion in the absence of an external mechanical force 
has been discussed in terms of the Marangoni effect in 
which a liquid droplet is driven by a surface tension gra- 
dient |13l 114] . The non-uniform surface tension can be 
controlled by an field variable such as temperature and a 
chemical (typically surfactants) concentration 15j. The 
mechanism is that the gradient induces convective flow 
inside and outside of a droplet, which leads to motion of 
the droplet itself. Similar flow and resulting motion are 
observed for a solid particle in phoretic phenomena such 
as thermophoresis [23 [T7]. In both systems, objects are 
swimming in a fluid. 

The velocity of the above-mentioned motion is reason- 
ably well described using linear theories [T31 [16j HI] • This 
implies that the direction of motion is determined by 
some asymmetry in the system such as a temperature 
gradient (and/or a concentration gradient). In the case 
of solid, an asymmetric particle has recently been created 
by coating half of its surface with a different material. 
Using this so-called Janus particle, the motion along a 
gradient created by the particle itself, which is referred 
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to as self-phoresis, was realized [19-21 . The asymmetric 
field in this case is not given externally but is created 
by consuming the energy supplied uniformly from out- 
side. Nevertheless, linear theory still works sufficiently 
well since the particle has inherently asymmetric surface 
properties. 

In contrast to the solid particle, fluid droplets are dy- 
namic and their surface properties cannot be fixed due to 
internal diffusion. Motion in an isotropic system cannot 
be described using a linear approach; it requires symme- 
try breaking arising from a nonlinear term |22) . In fact, 
spontaneous motion has been discussed using reaction- 
diffusion equations, which are nonlinear partial differen- 
tial equations, and is called as drift instability or drift 
bifurcation [23l425j . Despite this, there has been few at- 
tempts to consider the mechanics and hydrodynamics of 
spontaneous motion. 

In the present work, we derive amplitude equations 
showing drift bifurcation from a set of equations for con- 
centration fields taking hydrodynamics into considera- 
tion. All of the coefficients have clear physical mean- 
ings, and can in principle be measured. Our study is in- 
spired by earlier pioneering works on the motion of reac- 
tive droplets pjHi [2Z| ■ While these studies mainly focused 
on linear stability and response to an external force, our 
purpose is to derive equations containing nonlinear terms 
and obtain the characteristic velocity of a droplet. 



II. MODEL 

We consider an axisymmetric system containing a 
spherical droplet in a fluid which has an inner and/or 
outer surfactant concentration of c(r,9), and a velocity 
field of v(r, 0) = (v r (r, 9), ve(r, 9)) in the co- moving frame 
with the droplet [28] . Near the critical point of drift bi- 
furcation, the velocity of the droplet is slow so that v(r, 9) 
can be described by low-Reynolds hydrodynamics, that 
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FIG. 1. (Color Online) Schematic illustration of the system 
in this study. Surfactants dissolve in the outer fluid and some 
are adsorbed at the interface between the inner and the outer 
fluids. These surfactants reduce the surface tension of the 
droplet. (A) No droplet motion occurs for an isotropic distri- 
bution of surfactants. (B) When the surfactant distribution 
becomes asymmetric, the flow (thin red arrows) occurs and 
the droplet starts to move in the direction of the thick black 
arrow. The flux of surfactants are shown in broken arrows. 
The background gradation represents surfactant concentra- 
tion. 



is, the Stokes equation 



77V v = Vp, 



(1) 



with the incompressible condition V ■ v = 0. 77 is the vis- 
cosity of the inner or outer fluid and p is the pressure. We 
assume a linear relationship between the concentration of 
surfactants at the interface T(8) and surface tension 



7(0)=7o+7 c r(0) 



(2) 



using the surface tension 70 without surfactants. The 
surfactant concentration at the interface can be expanded 
using Legendre polynomials as 



r(e,t) 
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A n (t)P n (cos6). 



(3) 



Here we restrict our attention to non-deformable 
droplets. We consider only the n = and n — 1 modes, 
and neglect the higher modes. The solution of the Stokes 
equation for Marangoni flow for a given surface tension 
with an arbitrary distribution has been derived |14l I29j . 
It can be seen that the velocity of a droplet is propor- 
tional to the first mode as 



where 
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(4) 
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and the subscripts "i" and "o" denote the inner and outer 
fluid, respectively. 7 C is the strength of surface activity. 
Since the surface tension is typically smaller for higher 
concentrations of surfactants at the interface, 7 C is nega- 
tive and accordingly u\ > 0. u\ determines the strength 
of the chemomechanical coupling; the flow field is sensi- 
tive to the anisotropy when is large. Stronger cou- 
pling can be found for surfactants with higher surface 
activity. 



The concentration of molecules adsorbed at the inter- 
face is in balance with the bulk concentration field near 
the interface due to the adsorption-desorption equilib- 



rium as 



aT{6) = c(R,0), 



(6) 



where a is interpreted as the inverse of Henry's constant 
Kh for adsorption equilibrium and has the dimensions of 
inverse length |30j . For a low surfactant concentration 
at the interface, a is simply described as kd/k a where 
k a and kd are the adsorption rate from bulk and desorp- 
tion rate from surface, respectively. For this reason, a is 
not dimensionless but has the dimension of length. For 
surfactants with higher surface activity a can be small, 
for instance, a ~ 10 _1 m _1 [30 . The concentration of 
surfactants at the interface can be expressed as [28] 
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where the surface derivative is defined as V s = 
{\/R)d/d9 for a sphere. D (Di and D Q ) and D s are the 
bulk and surface diffusion constants, respectively. The 
surfactants are reactive; molecules dissolved in the bulk 
are adsorbed onto the interface, and after a characteristic 
time kJ 1 they lose their surfactant functionality, for in- 
stance by decomposing into a head and a tail (see FigJI]). 
We describe this by a linear reaction — k s T with a con- 
sumption rate n s . In this model, we implicitly assume ad- 
dition and removal of surfactants at the interface, which 
depend on the divergence of the two-dimensional velocity 
fields. 

We derive the amplitude equation near the onset of 
drift instability where A\ ~ e is small so that u ~ e. Our 
goal is to obtain the equation for the first mode 



m- 



dAi 
~dt~ 



gA 1 +eF 1 (A 1 ) + e 2 F 2 {A 1 ) 



(8) 



with coefficients m and g, and some functions J 7 !, J 7 ^, • • • ■ 
Taking Q into consideration, this is equivalent to a 
Landau-type equation for the droplet velocity 



fhdu/dt — gu + eP\{u) + ^^{u) + 



(9) 



The basic idea is to eliminate the velocity and bulk con- 
centration fields in order to obtain a closed form of the 
equations for A\. 

We hereafter focus only on the surfactant concentra- 
tion in the outer fluid and therefore drop the subscript 
"o" . The two fluids under consideration could, for exam- 
ple, be water and oil, and the surfactants preferentially 
dissolve in either one or the other. The bulk concentra- 
tion can be expressed using the Helmholtz equation with 
advection, 



dc 

dt 



VC: 
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(10) 
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The model takes into account the supply of surfactants to 
the bulk in order to maintain a constant concentration Coo 
far from the interface. The time scale is given by n. We 



expand ( 10 ) around the critical point of drift instability; 



the velocity of the droplet, u, or the Peclet number Ru/D 
is set as a small parameter e. We can solve this equation 
perturbatively as 



c(r, 9) = Coo + c (0) (r, 9) + ec w (r, 9) + 



(11) 



with the boundary conditions at infinity c(oo,0) = Cqq 
and at the interface (see (JgJ) ) . For the orders of e° and e 
, ( 10 ) is expressed as 



rfcW 



= DVV°> - kc<°>, (12) 
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The resulting c(r) is then substituted back into Q . Due 
to the boundary condition pi), the solution of c(r, 9) con- 
tains the individual modes A n and coupled modes A n A m . 
The nonlinear time evolution equations of A n are then 
obtained (see (16) and (17)). 



A. uniform distribution 

We assume that the relaxation of the bulk concentra- 
tion field is fast. The zeroth order solution of (12) is 
then 



c (°)(r,i) = Mo(f)- Coo ) 



kp(r/X) 
k (R/\y 



(14) 



where k n (x) is an nth-order modified spherical Bessel 
function of the second kind [3T] . The result is plotted in 
Fig{2| A) . A steep gradient can be observed in the typical 
length scale A = y/D/ k. The gradient is sustained by 
surface reaction characterized by n s in For R A, 
the surface concentration is given by 



A 



k s X/D 



(15) 



leading to a gap Coo — aAo between the concentration 
near the surface and at infinity. Since this gap is pro- 
portional to k s , the concentration gradient is driven by 
surface reactions. 



B. Amplitude equations 



A weakly nonlinear analysis up to the order of e 3 shows 

^=-Mo + ^(c»-aAo)+A^, (16) 
at AtC 
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FIG. 2. (Color Online) Distribution of bulk concentration 
field. (A) Isotropic distribution when A = 0.5, a — 0.01, and 
ui/D — 0. (B) Anisotropic distribution with Ui/D — 2.0. 
The blue (dark grey) line shows c(r, 6 — 0) (front) and the 
red (light grey) line shows c(r, 9 — n) (rear) . The uniform 
distribution of (A) is shown in (B) as a dashed line. 



where the coefficients are 
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(18) 
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with the coefficients A a b being dependent only on A and 
R. The explicit forms of A a t are shown in the Appendix 
(see ( |A49[ )-( [A53| )). The critical point of the drift bifur- 
cation occurs when the first term on the right-hand side 
of (17) changes its sign; 
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For u\ < u*, a stationary state is stable whereas it be- 
comes destabilized and the droplet moves for u\ > u*. 
For a <C n s \/D, the steady-state velocity of the droplet 
is given by 



Ui 



(22) 



where the characteristic velocity is uq — x/DkR/X for 
i?> A>0.01i?. 

The instability can be explained as follows. First, small 
fluctuations in the surfactant concentration at the inter- 
face give rise to a small A\, which induces convective 
flow around the droplet. The flow then distorts the bulk 
concentration field through the advection term. Above 
the critical point, the distortion overcomes the relaxation 
due to diffusion and amplifies the first mode A\ leading to 
further flow and motion of the droplet. In fact, Fig. [2jB) 
shows that the gradient in the bulk concentration at the 
front of the droplet (relative to the direction of motion) 
is steeper than that at the rear. This steeper gradient 
causes a larger flux from the bulk to the surface, and 
thus leads to an inhomogeneous surface concentration. 
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Above the critical point, the velocity increases with u\ 
as in FigjSJ In actual experiments, the size of a droplet 
may be the suitable parameter to vary. We find that 
there is an optimal droplet size for producing the highest 
velocity (Fig(3j3) . The two critical radii R\ ~ D s /cooUiA 
and i?2 — c ooU\\ 2 1 (Xk s + Da) arise from two stabilizing 
factors: surface diffusion and surface reaction. Both of 
them are balanced with the effect of advection. The size 
range for efficient self-propulsion increases with u\. The 
time evolution of the first mode below the critical point 
can be expressed as A\ ~ e~*/ Trolax , where the relaxation 
time is 



(101 
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(23) 



which diverges at u\ — u*. 

In the linear term of (17), A^/it* = A^Cqo — aAo), 
which corresponds to ( A25[ ) with (A51), destabilizes the 
stationary state. The physical origin of the destabiliza- 
tion is motion of the droplet. This can be seen in the first 



bracket in the velocity in radial direction ( A4 1 , which 
leads to the destabilization term. The first term in the 
bracket — uPi(cos6) corresponds to translational motion 
of the droplet in the co-moving frame while the second 
term u(R/r) 3 Pi (cos 9) arises from convective flow around 
the droplet. We investigated the contributions from both 
terms separately, and found that two terms have opposite 
effects; the first term (translational motion) destabilizes 
the stationary state while the latter (convection) stabi- 
lizes the instability. The instability is realized because 
the former always has stronger effect. 




FIG. 3. (Color Online) Bifurcation diagram for spontaneous 
motion. The bifurcation parameters are chosen to be u\ (A) 
and R (B). 
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(24) 
(25) 
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t = K s t, f = r/R, t = k s /k, l s = \J D s j n s R 2 , and I — 
•J 1 D I 'kR 2 . The velocity field is also non-dimensionalized 
as u = u\A\, where U\ = (Dcoo/kk s R 2 )ui, and A n — 
(Rk s / ' Dca^An. The boundary condition is rewritten as 
oT = 2(1) + 1 with a = (D / n s R)a. We choose u% to be 
a bifurcation parameter, which induces instability above 
a certain threshold, r is assigned a small value of 0.04. 
We estimate the critical point from the relaxation time 



using (23) 



We estimate the critical point from the relaxation time 
above the transition with ( 23 ) . Since the time evolution 



of A\ decays exponentially, we estimate the relaxation 
time by fitting the semi-log plot of A\ as a function of 
time. From the ir-intercept of the plot of relaxation time 
as a function of u\, we obtain the value of u\ at the 
critical point. The critical point weakly depends on the 
number of mesh points; for instance, for I — 0.2 and l s — 
1.0, our theory predicts — 1.03 while the numerical 
results show u\ = 1.09 for N = 100. As the mesh number 
is increased, the estimated critical point becomes closer 
to the predicted value u\ = 1.08 for N = 200 and u\ = 
1.04 for N = 400. Nevertheless, Fig. [4] shows that the 
normalized plot using numerically estimated values does 
not depend on the number of mesh points. We have 
mainly used N — 100 for saving computational time and 
for earning data points. 

The numerical results show the concentration distribu- 
tion around a droplet moving in the left direction |32j . It 
can be seen that the concentration distribution around 
the droplet is asymmetric. The droplet is stationary for 
small ui whereas it moves when u\ becomes larger. Note 
that the direction of motion is determined by an initially 
introduced small noise, and is therefore random. The ve- 
locity normalized by uq is plotted against the distance 
from the critical point in Fig. [3] Near the critical point 
the slope has a value of 0.5, which is comparable to the 
analytical result (22). A bifurcation is observed both 



III. NUMERICAL SIMULATIONS 



with and without the surface advection term in ([7]). The 
characteristic velocities deviate slightly from the analyt- 
ical results for some choices of parameters when the sur- 
face advection is included. This may be due to the effects 
of higher modes. Without the surface advection, all of 
the data points lie on the same curve irrespective of the 
parameters used. 



Numerical simulations are performed using spherical 
coordinates for an axisymmetric three-dimensional sys- 
tem. Both the radial and angular directions are dis- 
cretized with N + 1 mesh points. It is convenient to 
use the non-dimensionalized form of equations (|7| and 



IV. SUMMARY AND REMARKS 

In summary, we derive amplitude equations for drift 
instability of a droplet driven by Marangoni flow. The 
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FIG. 4. (Color Online) Normalized velocity of a droplet 
without (A) and with (B) surface advection in (24 1. The 
slope of the line is 0.5. 



Our model does not necessarily require the presence 
of surfactants. For instance, a uniformly heated droplet 
or a droplet with a source of chemicals can be tractable 
as the same manner with appropriate limits: a — > 1, 
D s — > 0, and n s » 1. In this situation, ^ is equivalent 
to a boundary condition for flux in the concentration field 
[D a n ■ Vc — Dili ■ Vc] r=fl — k s c(R) — 0, where the first 
and second terms represent the flux from outside and 
inside of the droplet, respectively. Here, the surface con- 
centration r independent of c(r) does not exist. However, 
it is convenient to introduce a virtual surface concentra- 
tion because velocity fields arc essentially created by the 
concentrations at the surface (see Q). It should also 
be stressed that similar results can be obtained using a 
phase-field model without explicitly considering a surface 

Although we focus on an outer fluid, generalization of 
the models to include the inner concentration is straight- 
forward. We may also consider production rather than 
consumption of surfactants at an interface, in which case 
spontaneous motion is realized for 7 C < 0. In fact, spon- 
taneous motion has been observed for complexes of sur- 
factants and ions that exhibit lower surface activity than 
the surfactants alone [TT] . 
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Appendix A: derivation of Eqs.( 18 |-( 20 1 



critical point and the droplet velocity arc calculated ana- 
lytically, and good agreement is found with the results of 
numerical calculations. Our system is out of equilibrium 
due to the reaction at the interface by which the supplied 
energy is consumed (see Q). This reaction maintains a 
concentration gradient in the radial direction. An addi- 
tional key factor is the nonlinear advection term in the 
bulk concentration field, which leads to coupling between 
modes and breaks the symmetry of the system. The con- 
centration gradient in the radial direction as well as the 
flux of surfactants onto the interface then becomes asym- 
metric. This leads to a surface tension gradient which 
results in motion. By contrast, surface advection is not 
essential for motility. Despite the linear nature of the 
velocity fields associated with the Stokes equation, we 
show that the addition of a nonlinear term in the con- 
centration field can lead to steady motion in an isotropic 
system. Further studies are required in order to clar- 
ify the kinds of nonlinear effects that are necessary for 
motility. 



In this appendix, we give a detailed derivation of the 
coefficients A ab in the amplitude equations (11) and (12). 
The dimensional analysis show that the coefficients have 
the dimension of length; they are functions of A and R. 
Introducing length and time scales, L and r, the param- 
eters are scaled as D ~ L 2 /t, a ~ 1/L, Ai ~ 1/L 2 , 
u = u\Ai ~ L/t, and u\ ~ L 3 /t . Then the coeffi- 
cients of amplitude equations are expressed as A02 ~ L , 
A 03 - L, An ~ 1/L, Ai 2 - L°, Ai 3 - L, and A u ~ L 2 . 
Using the coefficients, the steady velocity of a droplet is 
obtained from (12) as 
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Later, we will find A 12 ~ X/R and A14 ~ X 4 /R 2 which 
leads to the characteristic droplet velocity under a <C 
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k s \/D as 



U ~ Uq = 



DkR 
A ' 



(A2) 



In order to obtain the concrete form of the coeffi- 
cients, the Helmholtz equation with nonlinear advection 
is solved neglecting time derivative in (6), 
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v° • Vc^ 
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(All) 



V 2 c- ^(c- Coo ) = 



v - Vc 



(A3) 



for the order of e 2 , and so on. Note that although we 
focus only on the outer concentration field, the inner 
concentration field yields essentially the same equations. 



where the velocity field in the co-moving frame with the Hereafter, we drop the subscript o . 



droplet is given explicitly here as [21 [53] 
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P„(cos( 
(A4) 



The solution of the zeroth-order equation ( A9 ) satis- 
fying the boundary condition (3) is given in (9) using 
nth-order modified spherical Bessel function of the sec- 
ond kind k n (x) = y/lj {-KX)K, n+ i /% (x) where K n {x) is the 
nth-order modified Bessel function of second kind |31j . 
At the first order in the expansion, we will solve 
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Pi (cos ( 
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This equation is the form of the inhomogeneous 
Helmholtz equation: 
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which corresponds to ip = c 1 -- ^ an d / = using Z-th 
order expansion in (A10) and (All ). The inhomogeneous 
term is expanded as 

oo 

D , a , f( l \r,9)=J2fi l) (r)P n (cos9). (AH) 

^n(cos0), 

(A6) The general solution yields [34^ 

^(r^) = X)^n*n(r/A)P„(coflff)+ / G(r, r')/(r, 0)d 3 r, 

(A15) 

here the Green's function satisfies l3Tl 

sV, 

d9 (_„ 1 



dP n (cos9) 
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A 2 



G(r,r') = -S(r-r') 



(A16) 



where P„(cos#) is the nth-degree Legendre polynomial 
and 



For the Helmholtz equation in three dimensions, the 
Green's function is given as 



7c 



2(r?i + « ) 
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Near the critical point of drift bifurcation, the velocity of 
the droplet is small and accordingly the advection term 
is small. The solution is expanded perturb ative ly as c = 
Coo +c'°J +c^ 2 ) + • • • and at each order (A3 1 becomes 



47r|r - r'| 

The inhomogeneous term /W at the first order in ex- 



pansion (A12) is expressed as 

m e -(r-R)/x / ^3 

A (1) (r) = - „ 2 P(l- 



# o vV ) - KoC (°) = 
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——- Ai(coo - aA)) — , 
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and /„ = for n^l. The boundary condition (3) is which is shown in Fig. 2(B). Note that without the as- 
given as 

c {1) {R,9) = aAiP x {cos6). 
The solution would be 
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where 7li(r) is 



fci(r/A) / Mr'^r'/Xyidr' 
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+k(r/X) / A(r')fci(r7A)r' 2 dr 



(A21) 



i n (x) is the nth-order modified spherical Bessel function 
of first kind i n {x) — 'it / '(2x)X n+1 / 2 (x) using the nth- 
order modified Bessel function of first kind I n (x). This 
function has a simple form 

1 f°° 
KW(R) = -i 1 (R/\) h(r')h(r'/X)r' 2 dr'. (A22) 

A JR. 

The flux is expressed as 
<9c«(i?) aA x k[(r/\) 
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i[(R/X) k[(R/X) 
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For R > A (A22) becomes 

nW(R)^l^M Coo -aA )^ 
and we obtain the flux as 
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AxPi(cos0). 
(A25) 



The concrete form of TZ^(r) is 



su mptio n of i? > A the second term inside the bracket 
of (A25) is replaced by jr^^p (coo — aA )ui, which is 
always positive. This implies that this term destabilizes 
the stationary state irrespective of the value of A, that 
is, k. 

Calculation of the higher order terms is tedious but 
straightforward. The second order term in bulk concen- 
tration field satisfies 



1 



where 



/ (2) M) = 



uyAi 
D 



f - 



f (2 \r, 



(A28) 



P 3 \ dc« 



-Pi 



uiAx 



1 



Or 

R 3 \ I <9c« dPi 



(A29) 



D \~ ' 2r 3 J r 89 d9 
This is decomposed as 

/ (2) M) = /< 2) (r)P (cos 0) + /fV)^2 (cos 0) (A30) 
using 

Pi(cos(9)Pi(cos6i) = ip (cos<9) + ^P 2 (cos6<), (A31) 



dPi (cos 0) dPi (cos ( 



d$ 



d.e 



^P o (cos0)- ^P 2 (cos< 



(A32) 



Since we focus on the zeroth and first modes and the 
boundary condition is 



the general solution is expressed as 
J2)t-a\- ^(2)/ D ^ fc o(r/A) 



(A33) 



{ r,e)=-n^HR)^^+ll^(r), (A34) 



where 



fc« (r) = [2(r - P) 2 P(2r + R) 

+6r(r - R)RX - 3r(r - 2P)A 2 - 3rA 3 ] ( Coo - aA a ) 

(A26) 



ko(r/X) I ti 2 \r')i (r'/X)r' 2 dr' 

R 

r(2) 



+io(r/X) / tf>(r')k (r'/X)r a dr> 



(A35) 



and the concentration is 

A 1 Re-^ R '>/ x 
4r 3 (P + A) 



Similar to ( |A25 1, the flux is expressed as 



c«(r 



[ARar{r + A) 



dcW{R) 



D 



(-3r 2 P 2 - 3rP 2 A + 2r 3 (P + A) + R 3 (R + A)) 



dr 



i' (R/X) k' (R/X) 



x(coo - aA )] Pi (cose), 



(A27) 



io(P/A) k (R/X) 
271^1 (R) 
X 



-P (cos6»). 



P (cos6») 
(A36) 



with 



with 



P^(P) = 



u\A\ 



8D(R + X) 480£> 2 A 4 (P + A) 
x [A (4P 6 + 2P 5 A - P 3 A 3 + 3P 2 A 4 - 9PA 5 + 30A 6 ' 

-8e 2R/x R e (R + A)r[2P/A]l ( Coo - aA ). 

(A37) 

T[x] is the Gamma function. Note that the concentra- 
tion at this order is uniform since the coupling of two A\ 
modes results in Aq mode. For R 3> A, it is known that 
expansion does not converge |31) . Nevertheless, trunca- 
tion at finite terms in the series of expansion gives better 
approximation. 

The similar calculation is applied for the third-order 
equation: 



1 

A 2 



c^(r,9) = -f^(r,9), 



(3), 



(A38) 



where 



/ (3) (r,( 



u x A x ( R?\dc^ 



(A39) 



The solution is expressed as 



c< 3 )( 



where 



Pi (cos^), 

(A40) 



h(r/X) I f[ 3) (r')h(r' /X)r' 2 dr' 

R 



+*i(r/A) / tfHr^hir'/XyZdr' 



(A41) 



with 



/< 3 )(r,0) =/ 1 w (r)P 1 (cos^) + f^'(r)P 3 (cos9). (A42) 



r(3) 



(3)/ 



The flux is calculated as 



dc^(R) 



Or 



i[(R/X) k[(R/X) 
ii{R/X) ~ ki{R/X) 

2P.( 3 )(P) , n . 
^Pi(cos0) 



TZ^(R) 



Pi (cos ( 



(A43) 



p (3) (p) 

_ iu\A\ aX 4 (R - X){R 2 - 3PA + 3A 2 )(2P 2 + 6PA + 3A 2 ) 2 
80P 2 P 6 (P + A)(P 2 + 3PA + 3A 2 ) 



4AIB? 



240P 3 A 5 (P + A) 
where 

877 6139 A 



Cie -2R/x _c 2 r[2P/A] ( 



- aAo) 
(A44) 



Ci s 
and 

C 2 



19305 



94 



38610 R 

4 



26617 (X 
4290 VP 



+ 



541417 fX 



15444 



(A45) 



7016 
19305 

+ 220 



1754 R 80728 A 



19305 A 6435 P 



490814 / A 
6435 VP 



312 



(A46) 



We have used the integral including the Gamma function 

i r n+i A n+1 

r?r[0, n/X]dn = — r[o, r/A] + -— - I> + 1, r/A] 

n + 1 n + 1 

(A47) 

for n 7^ —1. In the limit of A — > 0, the second term 



of ( |A44| | becomes 1071A 6 w?Af /(128P 3 P). For the finite 
value of A, as mentioned above, the number of terms nec- 
essary for better approximation of the Gamma function 
depends on the value of A. For P ^> A > 0.01P, we have 



confirmed numerically ( A44 ) is well approximated by 
3u 2 AlaX 4 5ufAfX 5 



P (3) (P) 



20P 2 P 2 



6P 3 P 2 



(coo-aAo). (A48) 



The solution of c(r) is plugged into Ddc/dr in (4) and 
we obtain the set of amplitude equations (11) and (12). 
The coefficients are given as 



A 



Al2 



A 



A 



4(P + A)' 

03 - 32R2 , 

3PA 

4(P + A) 2 ' 

1 ' ~ 10P 2 ' 
5A 4 



14 



3P 2 ' 



(A49) 

(A50) 
(A51) 

(A52) 
(A53) 
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